%% calibrates the parameters of each region by matching to statistics from Figure 7

n_sim   = 100000;
Ty      = 20;
T       = ceil(Ty/dt);


%lower guarantees and higher leisure chi reduce the exit rate;
%in order to have a higher investment rate you need to increase r

if south==1
    sigma      = floor(100*sigma_south)/100;%when sigma is small the value function is very flat at small b and it has a hard time solving the problem; in particular the second derivative of v is zero
    varphi     = 0.0;
    delta      = delta_s; %if put delta=0, need to increase r
    coupon     = r+delta-varphi*delta;     %copupon
    A          = A_s;
    m_ss       = m_s_ss;
    y          = y_s;

    rho_vec      = rho;%this cannot be too small>=0.25 %calibrate to elasticity of investment to assets
    chi_vec     = nu/A*(0.01:0.01:0.05);%this has to be large enough so that there is a region with financial constraint and no default
    zeta_vec = 0.003:0.001:0.005;
    phi_vec      = 0.48:0.01:0.53;

    b_targ    = b_targ_south;
    exit_targ = bankrpcy_targ_south;
    exit_targ0= exit_targ_south;
    z_targ    = size_targ_south;
    str_region="South"
else
    sigma      = floor(100*sigma_north)/100;%when sigma is small the value function is very flat at small b and it has a hard time solving the problem; in particular the second derivative of v is zero
    varphi     = 0.0;
    delta      = delta_n; %if put delta=0, need to increase r
    coupon     = r+delta-varphi*delta;     %copupon
    A          = A_n;
    m_ss       = m_n_ss;
    y          = y_n;
    rho_vec      = rho;%this cannot be too small>=0.25 %calibrate to elasticity of investment to assets
    chi_vec     = nu/A*(0.02:0.01:0.06);%this has to be large enough so that there is a region with financial constraint and no default
    zeta_vec = [0.0075,0.01,0.0125];
    phi_vec      = 0.45:0.01:0.5;
    b_targ    = b_targ_north;
    exit_targ = bankrpcy_targ_north;
    exit_targ0= exit_targ_north;
    z_targ    = size_targ_north;
    str_region="North"
end

n_it = 0;
for i1=1:length(zeta_vec)
    for i2=1:length(chi_vec)
        for i3=1:length(phi_vec)
            for i4=1:length(rho_vec)
                n_it=n_it+1;
                phi_all(n_it)     = phi_vec(i3);
                chi_all(n_it)     = chi_vec(i2);
                zeta_all(n_it) = zeta_vec(i1);
                rho_all(n_it)      = rho_vec(i4);
            end
        end
    end
end
n_it

%% step 1
for ij=1:n_it;
    rho      = rho_all(ij);
    zeta = zeta_all(ij);
    chi     = chi_all(ij)*A/nu;
    phi      = phi_all(ij);
    a        = (a_targ-(1-phi)*varphi)/phi;

    inde_stop=0;
    get_firm_policy;
    %save policies
    pp_ell_loop(:,ij)  =pp_ell;
    pp_drift_loop(:,ij)=pp_drift;
    b_bar_loop(ij)     =b_bar;
    a_vec(ij)          =a;
end

calib=1;b_surv_mat=[];exit_prob_mat=[];z_surv_mat=[];
for ij=1:n_it
    ij
    rho      = rho_all(ij);
    zeta = zeta_all(ij);
    chi     = chi_all(ij)*A/nu;
    phi       = phi_all(ij);
    a          = a_vec(ij);

    pp_ell  =pp_ell_loop(:,ij);
    pp_drift=pp_drift_loop(:,ij);
    b_bar   =b_bar_loop(ij);

    b_up    = b_bar*0.75;
    b_md    = b_bar*0.25;
    q1_hig= (b_targ(1)-b_md)/(b_up-b_md);
    q_hig = q1_hig*b_up/b_md/(1-q1_hig*(1-b_up/b_md));
    q_mid = 1-q_hig;
    zmid    = 1/(1-q_hig+q_hig*b_md/b_up);
    zlow    = zmid*b_md/b_up;
    b_dn    = b_targ(1);zhig=1; q_low = 0.0;

    get_firm_simulations;

    b_surv_mat(:,ij)   = b_surv;
    exit_prob_mat(:,ij)= exit_surv;
    z_surv_mat(:,ij)   = z_surv;
end


str1="calib_";
str_save=append(path_save,str1,str_region,str_outcome,".mat");
filename=cell2mat(str_save);
save(filename, 'a_targ','exit_targ','b_targ','z_targ','chi_all','phi_all','zeta_all','rho_all','b_surv_mat' ,'exit_prob_mat','z_surv_mat','a_vec')

inde_ex=12:14;
inde_dz=12:14;
loss1 =   abs(squeeze(mean(exit_prob_mat(inde_ex,:)-exit_targ(inde_ex),1)))
loss2 =   abs(squeeze(mean(b_surv_mat(inde_ex,:)-b_targ(inde_ex),1)))
loss3 =   abs(squeeze(mean(z_surv_mat(inde_dz+1,:)./z_surv_mat(inde_dz,:)-z_targ(inde_dz+1)./z_targ(inde_dz),1)))
if south==1
    loss  =   loss1+loss2+loss3;
else
    loss  =   loss1+loss2+loss3;
end
[minM idx] = min(loss);

rho      = rho_all(idx);
zeta = zeta_all(idx);
chi     = chi_all(idx)*A/nu;
phi      = phi_all(idx);
a        = (a_targ-(1-phi)*varphi)/phi;
get_firm_policy;

str1="params_";
str2="_1";
str_save=append(path_save,str1,str_region,str2,str_outcome,".mat");
filename=cell2mat(str_save);
save(filename, 'delta', 'chi','zeta','sigma','varphi','rho','dr','b_bar','r','phi','a')

%% step 2: find the initial debt distribution
Ty      = 2;
T       = ceil(Ty/dt);
n_sim   = 50000;

share_targ_0 = b_targ(1)/b_bar;
if south==1
    share_targ_vec  = share_targ_0;
    bup_vec         = 0.2:0.005:0.3;
    bmd_vec         = -0.1:-0.005:-0.2;
else
    share_targ_vec  = share_targ_0;
    bup_vec         = 0.1:0.005:0.15;
    bmd_vec         = -0.01:-0.005:-0.1;
end

b_surv_mat_0   = zeros(Ty,length(bup_vec),length(bmd_vec));
exit_prob_mat_0= zeros(Ty,length(bup_vec),length(bmd_vec));
z_surv_mat_0   = zeros(Ty,length(bup_vec),length(bmd_vec));
b_0=zeros(length(bup_vec),length(bmd_vec));
for is=1:length(bmd_vec)
    share_targ=share_targ_0;
    b_targ_j=share_targ*b_bar;
    b_md    = (bmd_vec(is)+share_targ)*b_bar;
    for ij=1:length(bup_vec)
        ij
        tic
        inde_stop=0;
        b_up    = (bup_vec(ij)+share_targ)*b_bar;
        if b_up==b_md,q1_hig=0.5;
        else
            q1_hig= ( b_targ_j-b_md)/(b_up-b_md);
        end
        q_hig = q1_hig*b_up/(q1_hig*b_up+(1-q1_hig)*b_md);
        q_mid = 1-q_hig;
        zmid    = (1/b_md)/((1-q_hig)/b_md + q_hig/b_up);
        zlow    = zmid*b_md/b_up;

        zhig    = 1;
        q_low = 0;
        b_dn    =  b_targ_j;

        if zhig<0,warning('zhig'),inde_stop=1;end
        b_0(ij,is)=b_targ_j;
        if inde_stop==0 && b_up>=b_md && q_hig<.65 && q_hig>0.0

            get_firm_simulations;
            b_surv_mat_0(:,ij,is)   =b_surv;
            z_surv_mat_0(:,ij,is)   =z_surv;
            exit_prob_mat_0(:,ij,is)=exit_surv;
        else
            b_surv_mat_0(:,ij,is)   =1000*ones(Ty,1);
            z_surv_mat_0(:,ij,is)   =1000*ones(Ty,1);
            exit_prob_mat_0(:,ij,is)=1000*ones(Ty,1);
        end

    end
    toc
end

str1="calib_";
str2="_0";
str_save=append(path_save,str1,str_region,str2,str_outcome,".mat");
filename=cell2mat(str_save);
save(filename, 'b_surv_mat_0','z_surv_mat_0', 'exit_prob_mat_0', 'bup_vec', 'bmd_vec','b_0')

inde_val=1:2;
loss1=squeeze(max(abs(exit_prob_mat_0(inde_val,:,:)+delta-exit_targ0(inde_val)),[],1))+100*abs(b_0-b_targ(1));
loss2=squeeze(max(abs(exit_prob_mat_0(inde_val,:,:)-exit_targ(inde_val)),[],1))+100*abs(b_0-b_targ(1));
loss = loss1;
b_up_v    = b_bar*(bup_vec+share_targ_vec);%debt to value added (Z*A)
b_md_v    = b_bar*(bmd_vec+share_targ_vec);%debt to value added (Z*A)
[minM idx]    = min(loss(:));
[i1 i2]       = ind2sub(size(loss),idx);
% find values for province calibration of debt and productvity at entry
b_up    = b_bar*(bup_vec(i1)+share_targ_vec);%debt to value added (Z*A)
b_md    = b_bar*(bmd_vec(i2)+share_targ_vec);%debt to value added (Z*A)
b_targ_match=share_targ_vec*b_bar;
exit_prob_mat_0(inde_val,i1,i2)+delta- exit_targ0(inde_val)
exit_prob_mat_0(inde_val,i1,i2)- exit_targ(inde_val)

q1_hig= (b_targ_match-b_md)/(b_up-b_md);
q_hig = q1_hig*b_up/(q1_hig*b_up+(1-q1_hig)*b_md);
q_mid = 1-q_hig;
zmid    = (1/b_md)/((1-q_hig)/b_md + q_hig/b_up);
zlow    = zmid*b_md/b_up;

%find second chi of setting up the firm
v0_hig = ppval(pp_v,b_up);
v0_low = ppval(pp_v,b_dn);
v0_mid = ppval(pp_v,b_md);
x0_hig = ppval(pp_x,b_up);
x0_low = ppval(pp_x,b_dn);
x0_mid = ppval(pp_x,b_md);
Ev0    = A*(zlow*v0_hig*q_hig+zhig*v0_low*q_low+zmid*v0_mid*q_mid);%value of equity
B      = b_up*A*zlow;
vdebt_0= B*(x0_hig*q_hig+x0_low*q_low+x0_mid*q_mid);%cash flow from borrowing
str1="params_";
str2="_3";
str_save=append(path_save,str1,str_region,str2,str_outcome,".mat");
filename=cell2mat(str_save);
save(filename,  'b_up', 'b_dn', 'b_md', 'zhig', 'zmid', 'zlow','q_hig','q_mid','q_low', 'vdebt_0','Ev0','B')

%% step 3: simulate the entire age profile
Ty      = 20;
T       = ceil(Ty/dt);
n_sim   = 200000;

get_firm_simulations;
b_surv_mat_all   =[b_targ_match;b_surv(1:end-1)];
z_surv_mat_all   =z_surv;
exit_prob_mat_all=exit_surv;
str1="profile_";
str2="_all";
str_save=append(path_save,str1,str_region,str2,str_outcome,".mat");
filename=cell2mat(str_save);
save(filename,  'b_surv_mat_all', 'z_surv_mat_all', 'exit_prob_mat_all')

%% step 4: find the other parameters
dz=0.1;
z_max=15;
b_vec_mat_l = [-4.6:db:log(1+db),db];
z_vec       = log(max(0.05,dz)):dz:log(z_max);
db_mat      = db;
N_z =length(z_vec);
N_b =length(b_vec_mat_l);
Ntot=N_z*N_b;
str3="/codes/";
path_codes=append(root,str3);
addpath(path_codes)
get_matrices;

%% rescale fot the value added per province A_s
b_vec_mat   = exp(b_vec_mat_l+log(b_bar));
b_mat       = ones(length(z_vec),1)*(b_vec_mat);
bb_vec_mat  = reshape(b_mat,N_z*N_b,1);
ZZ          = reshape(exp(z_vec)'*ones(1,N_b),Ntot,1);
b           = reshape(ones(N_z,1)*(b_vec_mat),Ntot,1);

% find invariant distribution (b,z)
z_up=zhig;z_dn=zlow;z_md=zmid;inde_dilution=1;incid=1;policy_expe=0;
get_distribution_bZ;

avgz     = (ZZ'*f_dist)/sum(f_dist);
ret      = coupon./ppval(pp_x,min(b_bar,bb_vec_mat));
avgr     = (ret'*f_dist)/sum(f_dist)-delta

%rescale everything for firm value added and productivity (that's how data has been scaled)
z0       = y/(A*avgz*m_ss);%scaling factor of demand to match the level of output at the given A_s
zavg     = z0*avgz;
Z        = z0*ZZ;
B_0      = z0*B;
varzeta  = z0*Ev0*(r+delta);%this is in units of numeraire (final output)
k       = z0*vdebt_0+varzeta/(r+delta);% this is chi in units of numeraire (final output)
tau       = 0;

pars=[q_hig,zlow,q_low,zhig,q_mid,zmid,k,varzeta,rbar,tau,r,z0,delta,delta,A/nu-chi,chi,coupon,rho,sigma,phi,a,varphi,nu,1,y,m_ss,b_up,b_dn,b_md,B_0,zeta,theta,l_s,A,zavg];

str1="params_";
str2="_all";
str_save=append(path_save,str1,str_region,str2,str_outcome,".mat");
filename=cell2mat(str_save);
save(filename,  'pars','b_bar','pp_v','pp_x','pp_ell','pp_drift','pp_inv','Z')
